## R code for East Germany JPR paper Cox models Tables 1-5

load("~/Dropbox/Current projects/GDR/Spatial diffusion/Data/Modeling data 2.RData")

library(survival)
library(Hmisc)

## set stricter convergence criteria
coxph.control(eps = sqrt(.Machine$double.eps), toler.chol = .Machine$double.eps,
iter.max = 1000, toler.inf = sqrt(.Machine$double.eps), outer.max = 1000)


###################################
## Baseline models; Results Table 1
###################################
res1 <- matrix(NA,23,6)

m1a <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f + 

    		strata(eventnumbnumb),
    		data = protest6,
    		method = c("efron"),
    		model = TRUE)

summary(m1a)
res1[1:2,2] <- summary(m1a)$coef[1,c(1,3)]
res1[21,2] <- m1a$nevent
res1[22,2] <- m1a$n
res1[23,2] <- m1a$loglik[2]



m1b <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb),
    		data = protest6,
    		method = c("efron"),
    		model = TRUE)

summary(m1b)
res1[1:2,3] <- summary(m1b)$coef[1,c(1,3)]
res1[21,3] <- m1b$nevent
res1[22,3] <- m1b$n
res1[23,3] <- m1b$loglik[2]



m1c <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +
			
    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"),
    		model = TRUE)

summary(m1c)
res1[1:2,4] <- c(m1c$coef[1],sqrt(m1c$var[1,1]))
res1[21,4] <- m1c$nevent
res1[22,4] <- m1c$n
res1[23,4] <- m1c$history[[1]]$c.loglik



m1d <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +
			itm_dresden2f * wall + 

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m1d)
res1[1:2,5] <- c(m1d$coef[1],sqrt(m1d$var[1,1]))
res1[3:4,5] <- c(m1d$coef[2],sqrt(m1d$var[2,2]))
res1[5:6,5] <- c(m1d$coef[234],sqrt(m1d$var[234,234]))
res1[21,5] <- m1d$nevent
res1[22,5] <- m1d$n
res1[23,5] <- m1d$history[[1]]$c.loglik


## LR test
temp <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +
			wall + 

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

## LR test
LL.r <- temp$history[[1]]$c.loglik
LL.u <- m1d$history[[1]]$c.loglik
1 - pchisq((2*(LL.u - LL.r)), df = 1)



m1e <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +
			itm_dresden2f * thruoct18 + 
			itm_dresden2f * lateoct + 
			itm_dresden2f * nov1to9 + 

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m1e)
res1[1:2,6] <- c(m1e$coef[1],sqrt(m1e$var[1,1]))

res1[7:8,6] <- c(m1e$coef[2],sqrt(m1e$var[2,2]))
res1[9:10,6] <- c(m1e$coef[3],sqrt(m1e$var[3,3]))
res1[11:12,6] <- c(m1e$coef[4],sqrt(m1e$var[4,4]))

res1[13:14,6] <- c(m1e$coef[236],sqrt(m1e$var[236,236]))
res1[15:16,6] <- c(m1e$coef[237],sqrt(m1e$var[237,237]))
res1[17:18,6] <- c(m1e$coef[238],sqrt(m1e$var[238,238]))

res1[21,6] <- m1e$nevent
res1[22,6] <- m1e$n
res1[23,6] <- m1e$history[[1]]$c.loglik


res1a <- round(res1,3)
res1a[21:23,] <- round(res1a[21:23,])

latex(res1a, file = "")


## LR test
temp <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +
			thruoct18 + 
			lateoct + 
			nov1to9 + 

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))
    		
## LR test
LL.r <- temp$history[[1]]$c.loglik
LL.u <- m1e$history[[1]]$c.loglik
1 - pchisq((2*(LL.u - LL.r)), df = 3)



################
## results graph
################
out <- matrix(NA,3,9)
colnames(out) <- c("(1)","(2)","(3)","(4)Before","(4)After",
				"(5)Oct18","(5)LateOct","(5)EarlyNov","(5):End")
out[1,1:3] <- exp(res1[1,2:4])

out[2,1] <- exp(m1a$coef[1] - sqrt(m1a$var[1,1]) * qnorm(0.975))
out[3,1] <- exp(m1a$coef[1] + sqrt(m1a$var[1,1]) * qnorm(0.975))

out[2,2] <- exp(m1b$coef[1] - sqrt(m1b$var[1,1]) * qnorm(0.975))
out[3,2] <- exp(m1b$coef[1] + sqrt(m1b$var[1,1]) * qnorm(0.975))

out[2,3] <- exp(m1c$coef[1] - sqrt(m1c$var[1,1]) * qnorm(0.975))
out[3,3] <- exp(m1c$coef[1] + sqrt(m1c$var[1,1]) * qnorm(0.975))

# model 4
out[1,4] <- exp(m1d$coef[1] + m1d$coef[234])
out[1,5] <- exp(m1d$coef[1])

out[2,5] <- exp(m1d$coef[1] - sqrt(m1d$var[1,1]) * qnorm(0.975))
out[3,5] <- exp(m1d$coef[1] + sqrt(m1d$var[1,1]) * qnorm(0.975))

se <- sqrt(m1d$var[1,1] + m1d$var[234,234] + 2*m1d$var[1,234])
out[2,4] <- exp(m1d$coef[1] + m1d$coef[234] - se * qnorm(0.975))
out[3,4] <- exp(m1d$coef[1] + m1d$coef[234] + se * qnorm(0.975))

# model 5
out[1,6] <- exp(m1e$coef[1] + m1e$coef[236])
out[1,7] <- exp(m1e$coef[1] + m1e$coef[237])
out[1,8] <- exp(m1e$coef[1] + m1e$coef[238])
out[1,9] <- exp(m1e$coef[1])

out[2,9] <- exp(m1e$coef[1] - sqrt(m1e$var[1,1]) * qnorm(0.975))
out[3,9] <- exp(m1e$coef[1] + sqrt(m1e$var[1,1]) * qnorm(0.975))

se <- sqrt(m1e$var[1,1] + m1e$var[236,236] + 2*m1e$var[1,236])
out[2,6] <- exp(m1e$coef[1] + m1e$coef[236] - se * qnorm(0.975))
out[3,6] <- exp(m1e$coef[1] + m1e$coef[236] + se * qnorm(0.975))

se <- sqrt(m1e$var[1,1] + m1e$var[237,237] + 2*m1e$var[1,237])
out[2,7] <- exp(m1e$coef[1] + m1e$coef[237] - se * qnorm(0.975))
out[3,7] <- exp(m1e$coef[1] + m1e$coef[237] + se * qnorm(0.975))

se <- sqrt(m1e$var[1,1] + m1e$var[238,238] + 2*m1e$var[1,238])
out[2,8] <- exp(m1e$coef[1] + m1e$coef[238] - se * qnorm(0.975))
out[3,8] <- exp(m1e$coef[1] + m1e$coef[238] + se * qnorm(0.975))


plot(	seq(from = 1, to = 5.3, length.out = 6),
		seq(from = 0, to = 2, length.out = 6),
		type = "n", xaxt = "n",
		ylab = "WGTV effect", cex.lab = 1.3, xlab = "")

abline(h = 1, col = "grey", lwd = 2, lty = 3)

lines(x = c(1,1), y = c(out[1,1]-.03, out[2,1]), lwd = 4, col = "darkorange", lend = 1)
lines(x = c(1,1), y = c(out[1,1]+.03, out[3,1]), lwd = 4, col = "darkorange", lend = 1)
points(x = 1, y = out[1,1], pch = 20, col = "black", cex = 1)

lines(x = c(2,2), y = c(out[1,2]-.03, out[2,2]), lwd = 4, col = "darkorange", lend = 1)
lines(x = c(2,2), y = c(out[1,2]+.03, out[3,2]), lwd = 4, col = "darkorange", lend = 1)
points(x = 2, y = out[1,2], pch = 20, col = "black", cex = 1)

lines(x = c(3,3), y = c(out[1,3]-.03, out[2,3]), lwd = 4, col = "darkorange", lend = 1)
lines(x = c(3,3), y = c(out[1,3]+.03, out[3,3]), lwd = 4, col = "darkorange", lend = 1)
points(x = 3, y = out[1,3], pch = 20, col = "black", cex = 1)

# model 4
lines(x = c(3.9,3.9), y = c(out[1,4]-.03, out[2,4]), lwd = 4, col = "darkorange",
	lend = 1)
lines(x = c(3.9,3.9), y = c(out[1,4]+.03, out[3,4]), lwd = 4, col = "darkorange",
	lend = 1)
points(x = 3.9, y = out[1,4], pch = 20, col = "black", cex = 1)

lines(x = c(4.1,4.1), y = c(out[1,5]-.03, out[2,5]), lwd = 4, col = "darkorange",
	lend = 1)
lines(x = c(4.1,4.1), y = c(out[1,5]+.03, out[3,5]), lwd = 4, col = "darkorange",
	lend = 1)
points(x = 4.1, y = out[1,5], pch = 20, col = "black", cex = 1)


# model 5
lines(x = c(4.7,4.7), y = c(out[1,6]-.03, out[2,6]), lwd = 4, col = "darkorange",
	lend = 1)
lines(x = c(4.7,4.7), y = c(out[1,6]+.03, out[3,6]), lwd = 4, col = "darkorange",
	lend = 1)
points(x = 4.7, y = out[1,6], pch = 20, col = "black", cex = 1)

lines(x = c(4.9,4.9), y = c(out[1,7]-.03, out[2,7]), lwd = 4, col = "darkorange",
	lend = 1)
lines(x = c(4.9,4.9), y = c(out[1,7]+.03, out[3,7]), lwd = 4, col = "darkorange",
	lend = 1)
points(x = 4.9, y = out[1,7], pch = 20, col = "black", cex = 1)

lines(x = c(5.1,5.1), y = c(out[1,8]-.03, out[2,8]), lwd = 4, col = "darkorange",
	lend = 1)
lines(x = c(5.1,5.1), y = c(out[1,8]+.03, out[3,8]), lwd = 4, col = "darkorange",
	lend = 1)
points(x = 5.1, y = out[1,8], pch = 20, col = "black", cex = 1)

lines(x = c(5.3,5.3), y = c(out[1,9]-.03, out[2,9]), lwd = 4, col = "darkorange",
	lend = 1)
lines(x = c(5.3,5.3), y = c(out[1,9]+.03, out[3,9]), lwd = 4, col = "darkorange",
	lend = 1)
points(x = 5.3, y = out[1,9], pch = 20, col = "black", cex = 1)


axis(side = 1, at = 1:5, labels = c("(1)","(2)","(3)","(4)", "(5)"), cex.axis = 1,
		tick = FALSE)

text(x = 3.9, y = 0.3, "Before fall of Berlin Wall", srt = 90, cex = .75)
text(x = 4.1, y = 0.3, "After fall of Berlin Wall", srt = 90, cex = .75)

text(x = 4.7, y = 0.3, "Throughout October 18", srt = 90, cex = .75)
text(x = 4.9, y = 0.3, "Late October", srt = 90, cex = .75)
text(x = 5.1, y = 0.3, "Early November", srt = 90, cex = .75)
text(x = 5.3, y = 0.3, "After fall of Berlin Wall", srt = 90, cex = .75)



#####################
## Robustness Table I
#####################
res2 <- matrix(NA,13,6)

m2a <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

			neighbor3 +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m2a)
res2[1:2,2] <- c(m2a$coef[1],sqrt(m2a$var[1,1]))
res2[11,2] <- m2a$nevent
res2[12,2] <- m2a$n
res2[13,2] <- m2a$history[[1]]$c.loglik



m2b <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

			neighbor7  +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m2b)
res2[1:2,3] <- c(m2b$coef[1],sqrt(m2b$var[1,1]))
res2[11,3] <- m2b$nevent
res2[12,3] <- m2b$n
res2[13,3] <- m2b$history[[1]]$c.loglik



m2c <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

			neighbor14  +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m2c)
res2[1:2,4] <- c(m2c$coef[1],sqrt(m2c$var[1,1]))
res2[11,4] <- m2c$nevent
res2[12,4] <- m2c$n
res2[13,4] <- m2c$history[[1]]$c.loglik



m2d <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

			neighbor3 + neighbor7 + neighbor14 +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m2d)
res2[1:2,5] <- c(m2d$coef[1],sqrt(m2d$var[1,1]))
res2[11,5] <- m2d$nevent
res2[12,5] <- m2d$n
res2[13,5] <- m2d$history[[1]]$c.loglik



m2e <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

			poland + trains +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m2e)
res2[1:2,6] <- c(m2e$coef[1],sqrt(m2e$var[1,1]))
res2[11,6] <- m2e$nevent
res2[12,6] <- m2e$n
res2[13,6] <- m2e$history[[1]]$c.loglik


res2a <- round(res2,3)
res2a[11:13,] <- round(res2a[11:13,])

latex(res2a, file = "")



################
## ROBUSTNESS II
################
res3 <- matrix(NA,21,9)

m3a <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_852f  +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m3a)
res3[1:2,2] <- c(m3a$coef[1],sqrt(m3a$var[1,1]))
res3[19,2] <- m3a$nevent
res3[20,2] <- m3a$n
res3[21,2] <- m3a$history[[1]]$c.loglik



m3b <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_825f  +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m3b)
res3[3:4,3] <- c(m3b$coef[1],sqrt(m3b$var[1,1]))
res3[19,3] <- m3b$nevent
res3[20,3] <- m3b$n
res3[21,3] <- m3b$history[[1]]$c.loglik



m3c <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_80f  +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m3c)
res3[5:6,4] <- c(m3c$coef[1],sqrt(m3c$var[1,1]))
res3[19,4] <- m3c$nevent
res3[20,4] <- m3c$n
res3[21,4] <- m3c$history[[1]]$c.loglik



m3d <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_prop_dresdenf   +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m3d)
res3[7:8,5] <- c(m3d$coef[1],sqrt(m3d$var[1,1]))
res3[19,5] <- m3d$nevent
res3[20,5] <- m3d$n
res3[21,5] <- m3d$history[[1]]$c.loglik



m3e <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_prop_85f   +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m3e)
res3[9:10,6] <- c(m3e$coef[1],sqrt(m3e$var[1,1]))
res3[19,6] <- m3e$nevent
res3[20,6] <- m3e$n
res3[21,6] <- m3e$history[[1]]$c.loglik



m3f <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_prop_825f   +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m3f)
res3[11:12,7] <- c(m3f$coef[1],sqrt(m3f$var[1,1]))
res3[19,7] <- m3f$nevent
res3[20,7] <- m3f$n
res3[21,7] <- m3f$history[[1]]$c.loglik



m3g <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_prop_80f   +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m3g)
res3[13:14,8] <- c(m3g$coef[1],sqrt(m3g$var[1,1]))
res3[19,8] <- m3g$nevent
res3[20,8] <- m3g$n
res3[21,8] <- m3g$history[[1]]$c.loglik



## recode tv5 to make it identical to measure used in Kern (2011)
protest6$tv5f[protest6$county == "Doebeln"] <- 1
protest6$tv5f[protest6$county == "Grimma"] <- 1
protest6$tv5f[protest6$county == "Freiberg"] <- 1

table(protest6$county,protest6$tv5f)
# note: One county is missing since it did not experience any protests

m3h <- coxph(formula = Surv(istart, ifinish, protest) ~
    		tv5f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		method = c("efron"))

summary(m3h)
res3[15:16,9] <- c(m3h$coef[1],sqrt(m3h$var[1,1]))
res3[19,9] <- m3h$nevent
res3[20,9] <- m3h$n
res3[21,9] <- m3h$history[[1]]$c.loglik


res3a <- round(res3,3)
res3a[19:21,] <- round(res3a[19:21,])

latex(res3a, file = "")



##################
### ROBUSTNESS III
##################
res4 <- matrix(NA,7,7)

## omit Berlin
sel<- protest6$county != "Berlin (Ost)"
table(sel)

m4a <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		subset = sel,
    		method = c("efron"))

summary(m4a)
res4[1:2,2] <- c(m4a$coef[1],sqrt(m4a$var[1,1]))
res4[5,2] <- m4a$nevent
res4[6,2] <- m4a$n
res4[7,2] <- m4a$history[[1]]$c.loglik



## only use counties in the Northeast as control group
sel<- 	protest6$county == "Bautzen" |
 		protest6$county == "Dresden (Stadt)" |   
		protest6$county == "Forst" |
		protest6$county == "Goerlitz (Stadt)" |
		protest6$county == "Goerlitz-Land" |
		protest6$county == "Hoyerswerda" |
		protest6$county == "Loebau" |
		protest6$county == "Niesky" |
		protest6$county == "Pirna" |
		protest6$county == "Weisswasser" |
		protest6$county == "Zittau"

sel <- !sel
table(sel)

m4b <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		subset = sel,
    		method = c("efron"))

summary(m4b)
res4[1:2,3] <- c(m4b$coef[1],sqrt(m4b$var[1,1]))
res4[5,3] <- m4b$nevent
res4[6,3] <- m4b$n
res4[7,3] <- m4b$history[[1]]$c.loglik



## only use counties in the Southeast as control group
sel2 <- protest6$itm_dresden2f == FALSE & !sel == FALSE
sel2 <- !sel2
table(sel2)
table(protest6$county[!sel2])

m4c <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		subset = sel2,
    		method = c("efron"))

summary(m4c)
res4[1:2,4] <- c(m4c$coef[1],sqrt(m4c$var[1,1]))
res4[5,4] <- m4c$nevent
res4[6,4] <- m4c$n
res4[7,4] <- m4c$history[[1]]$c.loglik

rm(sel,sel2)



## drop WGTV counties that neighbor non-WGTV counties
load("~/Dropbox/Current projects/GDR/Dave analysis/Omit counties RR.RData")

sel <- is.element(protest6$county, omit.counties)
table(sel)
table(protest6$county[sel])
sel <- !sel

m4d <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f +
    		exits_pop +
			protests53_cities +

    		strata(eventnumbnumb) +
    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		subset = sel,
    		method = c("efron"))

summary(m4d)
res4[1:2,5] <- c(m4d$coef[1],sqrt(m4d$var[1,1]))
res4[5,5] <- m4d$nevent
res4[6,5] <- m4d$n
res4[7,5] <- m4d$history[[1]]$c.loglik



## (e) copy and paste from balance ALL results.



## (f) only use time to first event
sel<- protest6$eventnumbnumb > 1 | (protest6$eventnumbnumb == 1 & protest6$protest == 0) 
sel <- !sel
table(sel)

m4f <- coxph(formula = Surv(istart, ifinish, protest) ~
    		itm_dresden2f +

			logpopulation_1988f + 
    		femalef100 +
    		collegef100 +
    		industryf100 +
    		working_agef100 +
    		density_1988f +
    		pop_changef +
    		agriculturef +
    		craftsconstructionf +
    		servicestransportf +
    		unskilledf +
    		skilledf +
    		housing_spacef +
    		bathroomf +
    		interior_toiletf +
    		modern_heatingf +
    		capita_per_mdf + 
    		capita_per_dentistf +
    		nitrogen_oxidesf +
    		sulfur_dioxidef +
    		respirable_dustf +
    		districtcapital2f + 
    		exits_pop +
			protests53_cities +

    		frailty(id, distribution = "gamma", method = "em", sparse = FALSE),
    		data = protest6,
    		subset = sel,
    		method = c("efron"))

summary(m4f)
res4[1:2,7] <- c(m4f$coef[1],sqrt(m4f$var[1,1]))
res4[5,7] <- m4f$nevent
res4[6,7] <- m4f$n
res4[7,7] <- m4f$history[[1]]$c.loglik


res4a <- round(res4,3)
res4a[5:7,] <- round(res4a[5:7,])

latex(res4a, file = "")













.


